Structure and correlations in two-dimensional 
Coulomb confined classical artificial atoms 
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The ordering of TV equally charged particles (— e) moving in two dimensions and confined by a 
Coulomb potential, resulting from a displaced positive charge Ze is discussed. This is a classical 
model system for atoms. We obtain the configurations of the charged particles which, depending on 
the value of TV and Z, may result in ring structures, hexagonal-type configurations and for N/Z « 1 
in an inner structure of particles which is separated by an outer ring of particles. For N/Z << 1 
the Hamiltonian of the parabolic confinement case is recovered. For N/Z ps 1 the configurations are 
very different from those found in the case of a parabolic confinement potential. A hydrodynamic 
analysis is presented in order to highlight the correlations effects. 

PACS numbers: 71.15 Pd, 73.21 La 



I. INTRODUCTION 



Quantum dots, or artificial atoms, have been a subject 
of intense theoretical and experimental studies in recent 
years jjj due to the occurrence of numerous interesting 
effects caused by electron correlation, such as e. g. Wigner 
crystallization, overcharging, and nontrivial behavior in 
a magnetic field. These electron correlations appear most 
clearly when the electron interaction and the confinement 
potential dominates over the kinetic energy of the system. 
This can be realized in e.g. quantum dots which have 
much smaller electron density as compared to real atoms. 

Correlation effects show up in an even more pro- 
nounced way in classical systems where the kinetic en- 
ergy is zero in the absence of thermal fluctuations. Two- 
dimensional (2D) classical dots confined by a parabolic 
potential were studied earlier and a table of Mendeleyev 
for such " artificial" atoms was constructed m . The clas- 
sical system which is more closely related to real atoms 
was studied in Ref. H where, as function of the strength 
of the confinement potential, surprising rich physics was 
observed like structural transitions, spontaneous symme- 
try breaking, and unbinding of particles, which is absent 
in parabolic confined dots. The confinement potential 
was of the Coulomb type, but in order to prevent the 
collapse of all electrons onto the "nucleus" (which we 
call 'impurity' in this case) the positive charge (Ze) was 
displaced a certain distance from the 2D plane where the 
electrons are moving in (see Fig. III). Note that our sys- 
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tern is related to the "superatom" system introduced by 
Watanabe and Inoshita Q . The superatom is a spherical 
modulation-doped heterojunction. In particular it is a 
quasi-atomic system which consists of a spherical donor- 
doped core and a surrounding impurity-free matrix with 
a larger electron affinity. The quantum mechanical elec- 
tron structure of this system was studied in S and it 
was found that due to the absence of the 1 /r singularity 
in the potential the ordering of the energy levels is dom- 
inated by the no-radial- node states, in contrast to real 
atoms where s and p states are dominant. 

In the present paper we present a systematic study 
of the system of Ref. || as a function of the number 
of particles (TV). In contrast to Ref. || we will not 
vary the strength of the confinement potential, as this 
was already presented in our earlier work || where we 
limited the numerical results to the case Z = N = 4, 
but we will vary Z and N. Furthermore, we will compare 
our numerical results with the results of a hydrodynamic, 
i.e. continuum, approach. In the latter case the electron 
density is taken as a fluid, i.e. there are no charge quanta. 
This is in contrast to our numerical simulation where 
the electrons are point particles with a fixed quantized 
charge value. The fact that charge is now distributed in 
packages of q = — e will introduce important correlation 
effects which is the central theme of the present work. 

Note also that the present classical study can serve 
as a zeroth order approach for more demanding quan- 
tum mechanical calculations. Recently, such classical 
calculations were used || as a starting point in order 
to construct better quantum wavefunctions, at least in 
the strong magnetic field limit. 

Besides the above mentioned analogy with real and 
artificial atoms which are inherently quantum mechani- 
cal there exist other experimental realized systems which 
behave purely classical and for which our study is rele- 
vant. Examples are charged colloidal suspensions where 
it was found recently that correlation effects between the 
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FIG. 1: Schematic view of the system. 



counter-ions can result into an overscreening and attrac- 
tion between like charged colloids Q. Our system is a 
simplified 2D model for those colloidal systems. 

The system under study can be realized experimen- 
tally in the system of electrons above liquid helium || by 
putting a positive localized charge in the substrate which 
supports the liquid helium. The equivalent quantum me- 
chanical system can be realized using low dimensional 
semiconductor structures with impurities, also called re- 
mote impurities, which are displaced a distance from a 
quantum well Q . In both cases it will be rather difficult 
to increase the number of positive charge Z beyond a few 
units. 

Another possible realization of our system is by bring- 
ing an AFM (atomic force microscope) tip close to a 2D 
electron gas. When this tip is charged positively it will 
induce a confinement potential very similar to the one 
studied in the present paper. The advantage of this ap- 
proach is that the charge Ze on the tip can be varied 
continuously by increasing the voltage on the tip. 

This paper is organized as follows. In section II we 
describe the mathematical model and our numerical ap- 
proach to obtain the configurations. The results of our 
numerical simulations are given in section III. A hydro- 
dynamic approach which neglects the correlation effects 
is presented in sections IV. The numerical simulation re- 
sults are compared with the results of the hydrodynamic 
approximation in section V which clearly brings about 
the importance of the charge correlation. Our conclu- 
sions are presented in section VI. 



II. THE MODEL 

We study a system with N negatively charged particles 
(— e), which we call here and further electrons, interacting 
through a repulsive Coulomb potential and moving in 
the sy-plane. The particles are kept together through a 
fixed positive charge (Ze) located at a distance a from 
the plane the particles are moving in (see Fig. 0). The 
total energy of this system is given by the Hamiltonian 
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Here the symbol e stands for the dielectric constant, and 
r = {x, y} is the two component position vector of the 2D 



(two dimensional) electron. For convenience, we express 
the electron energy in units of Eq = e 2 /ea and all the 
distances in units of a. This allows us to rewrite Eq. (|l|) 
in the following dimcnsionless form: 
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The ground state configurations of the two- 
dimensional system were obtained using the standard 
Metropolis algorithm jn^. The electrons are initially 
put in random positions within some circle and allowed 
to reach a steady state configuration after a number 
of simulation steps in the order of 10 5 . To check if 
the obtained configuration is stable, we calculated the 
frequencies of the normal modes of the system using 
the Householder diagonalization technique jl2). The 
configuration was taken as final when all frequencies of 
the normal modes are positive and the energy did not 
decrease further. The meta-stable states were avoided 
introducing a small temperature T = 10~ 4 , which 
was negligible and does not influence the simulation 
accuracy. 



III. STABLE CONFIGURATIONS 

As an example of our results we present in Fig. || the 
radial distribution of the electrons for a fixed positive 
charge Z — 50 for the impurity as function of the num- 
bers of electrons N in the dot. Clearly, two different 
electron distribution types can be distinguished. Namely, 
for small number of electrons N a ring structure arrange- 
ment for the electrons is observed which is similar to the 
one for a parabolic dot (compare with the configurations 
given in ||), for large N the outer electrons can form a 
ring which is clearly separated from the other electrons in 
the dot. This is more clearly seen in Fig. ^| where exam- 
ples of two different configurations (for small and large 
N) are presented. Note that in the case of large N the 
core electrons are arranged in nearly a triangular lattice 
which is characteristic for an infinite electron system. 

The very different type of configurations for small N 
and large N is a consequence of the fact that the screen- 
ing of the Coulomb center is a much more complicated 
problem as compared with the parabolic dot case. Now 
its behavior is controlled by two parameters, namely, the 
number of electrons N, which characterizes the discrete- 
ness of the charges taking participation in the screening, 
and the positive charge Z, which represents the strength 
of the confinement Coulomb potential. 

In the case of a relatively small number of electrons 
(n = N/Z <C 1) the confinement potential is strong, and 
the electrons are located close to the origin where the 
confinement potential can be replaced by the following 
approximate one: 
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FIG. 2: The radial position of the particles, n, as a function 
of the number of electrons in the dot, N, for a system with 
confinement positive charge Z = 50. The number of electrons 
is varied from N = 2 up to N = 50. 
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FIG. 3: Configurations for two values of the number of parti- 
cles N for a fixed Coulomb center with Z = 50. 



Now substituting it into Hamiltonian (|2j) , and scaling the 
variables 

H^{Z/2) 1/3 H, r ^ (2/Z) 1/3 r, (4) 
one obtains the Hamiltonian 
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FIG. 4: The energy per particle (a) and the maximum radius 
(6) of the Coulomb dot (solid circles) and the parabolic dot 
(open triangles) as a function of N. The confinement charge 
of the Coulomb dot is Z = 200. 



which, up to a non-essential energy shift, coincides with 
the Hamiltonian of a parabolic dot as considered in ref- 
erence In the above Z = 50 case we have the same 
parabolic dot like configurations up to N = 9, while for 
N > 10 new configurations as (3,7), (4,7) and so on ap- 
pear. For larger Z the similarity between parabolic dot 
and Coulomb center screening persists up to larger N 
values. 

The equivalence of Coulomb screening and the results 
for a parabolic dot at small N values is confirmed in 
Fig. ||, where the energy per particle and the maximum 
radius of a Coulomb dot confined by a positive impu- 
rity charge Z = 200 are shown as a function of N and 
compared with the results obtained for the parabolic dot 
^ More over, if one plots the differences of energy 
per particle (scaled by Z) as a function of the number 
of particles N, as is shown in Fig. [jj, we can see small 
cusps related to the so called "magic numbers" which are 
known to be an important feature of the configurations 
of the parabolic dots |T^| . 

In the opposite case of large electron numbers (n < 1), 
where the system is nearly neutral, one can expect to 
find features which are specific for our Coulomb screen- 
ing case and which are not found with a parabolic dot. 
In this asymptotic case, the Coulomb dot presents a low 
density system H which is under the influence of a non- 
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FIG. 5: Difference in the energy per particle per Z as a func- 
tion of N for a Coulomb dot with Z = 100. Examples of 
configurations where kinks are found in the energy curve are 
shown in the insets. 



homogeneous electric field. We already know that cor- 
relation effects in such a system is of most importance. 
That is why it is worth to take as a reference the hy- 
drodynamic approach, which is some mean field theory 
where no correlation effects are included. 



IV. HYDRODYNAMIC APPROACH 

In the hydrodynamic approach the electrons are de- 
scribed by the 2D density p(r). They create a potential 
4>{r, z) which obey the Poisson equation 



(r, z) — 4atp(r)5(z), 



(6) 



in the whole 3D space. Using the dimensionless variables 
introduced in Sec. H the electron density is measured 
in a~ 2 units and the potential in e/a units. Note we 
took into account the cylinder symmetry of the problem 
caused by the confining potential V(r) — —Zj\]r 2 + 1 
(see Eq. (||)). The Poisson equation can be replaced by 
the Laplace equation 



r,z)=0 
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everywhere outside the xy-plane together with the 
boundary condition 



dc/)(r, z) 



dz 



= 2np(r) 
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on this plane. 

The mathematical model is based on the statement 
that the electrons (located in the circular disk of ra- 
dius R) have the same energy in every point of the dot, 
namely, 



{V(r) - 0(r)} 



r<R 



p = const. 



(9) 



Here the symbol (f>(r) = <fi(r, 0) stands for the potential 
created by the 2D electrons in the disk, and the symbol 
p is the chemical potential of this electron system. 

The above equations have to be supplemented with one 
more expression, namely, the conservation of the total 
number of electrons 



N = 2tt 



rdrp{r). 



(10) 



The hydrodynamic approach is a kind of mean field 
theory where no correlation effects are included. For in- 
stance, the phenomena of overcharging Q cannot take 
place, and the condition n = N/Z < 1 is always satis- 
fied, namely, the number of electrons N attracted by the 
Coulomb center never exceeds its charge number Z. The 
specific case n = 1 is referred to as complete screening. 
In this case the analytical solution of the above equations 
can be obtained using the mirror charge technique, and 
it leads to the following result 
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with n cs — 1 and p cs = 0. This simple limiting case 
result is useful as a reference. 

In the general N < Z case we solved the hydrodynamic 
equations using the oblate spherical coordinates (0 < r < 
1, < <t < oo) which are defined by 



x = i?v/(cr 2 + l)(l-T 2 )cos6>, 
y = R^{a 2 + 1)(1 - T 2 )sin6>, 
z = Ro~t 



(12a) 
(12b) 
(12c) 



as it was done in Refs. |Tl], |l3| for the case of a parabolic 
dot with radius R. The solution of the Laplace equation 
(0) can be presented as an expansion 



(r, z) = $(r,cr) 



(13) 



in terms of the first Pi n and second Q2n kind Legendre 
polynomials. Thus, the potential created by the electrons 
on the disk (cr = 0, r < 1) can be presented as 
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Now taking int o ac count that on the disk we have ac- 
cording to Eq. ( po] ) 
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and satisfying the boundary condition (||) one gets the 
analogous expansion for the electron density 



p[r) = X (r) - C nLnP2n{T 



(16) 
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and r(a;) is the Gamma function. 

Next, we expand the Coulomb center potential on the 
disk into a series of Legendre polynomials as well 



V(r) = - 



Ry/l + 1/R 2 - T 2 
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and inserting it together with Eq. (|l4| ) into Eq. we 
obtain the final expression 
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for the electron density expansion coefficients C n . 

In order to define the chemical potential fj, we have to 
remember that the electron density has to be equal to 
zero at the free electron system boundary r = R. Thus, 
the electron density should satisfy the following condition 



2irR lim t X (t) = - V C n L n P 2n (0) 

T—>0 Z 4 
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FIG. 6: Hydro dynamic solution (solid curves) and its asymp- 
totic (dashed curves) for large R: 1 - number of particles 
(Eq. @) , 2 - chemical potential (Eq. @). 
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what finally enables us to define the chemical potential 
Z 



RL{ 



n=0 



n— 



and write down the following electron density expression 



X(r) = 



Z 



2irR 2 ', 



{P2n(r) - P 2 „(0)} . (22) 



Inserting the density expression (|l(]) into Eq. ( |io| ) one 
obtains the number of electrons 



N 2 ( A R 

n = ^ = n{ Ao+ z fl 



(23) 



taking part in the screening of the Coulomb center. 

Equations ([TJ|), (|TJ), (||), and (||) actually are the 
solution of the problem if the coefficients A n are known. 
These coefficients are calculated in Appendix |A| 

The density ( p2| ) enables us to calculate all properties 
of the electron system. For instance, the square of mean 
electron radius can be estimated as 

< r *>=i rw^* (|-i^). (24, 



FIG. 7: Hydrodynamic solution for the electron density for 
various dot radii. 



Due to the linear dependence of the confinement poten- 
tial ( |l8| ) , the chemical potential ( pl| ) and the density J22] ) 
on the charge number Z , it can be removed by a scaling, 
and the hydrodynamic problem is actually controlled by 
a single parameter (say, n or R) in contrast to the sys- 
tem of discrete electrons for which the two parameters (Z 
and N) were essential. Thus, the general solution can be 
presented by both chemical potential and density curves, 
as is shown in Fig. |s| by solid curves. The correspond- 
ing dashed curves indicate the asymptotic of the above 
quantities 



= 1 - 



Z 

R 2 ' 
4 

7ri?' 



(25a) 
(25b) 



which holds for a nearly neutral (n — > 1) electron sys- 
tem. In Fig. |?] the electron density for various dot radii 
is shown. We see that in the case of small radius (or 
equivalently large Z, i.e. n << 1) the electron den- 
sity becomes similar to the density in a parabolic dot 
(p ~ yfl — r 2 / R 2 ), while in the opposite large radius 
case it tends to the limiting density ( pd| ) for the com- 
pletely screened case. 
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FIG. 8: The scaled energy per particle (a) and chemical po- 
tential (b) obtained from our numerical simulations for dif- 
ferent values of the confinement charge Z are compared with 
the hydrodynamic result. 
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FIG. 9: The difference in energy (a) and chemical potential 
(b) between the hydrodynamic result and the simulation re- 
sult for Z = 100. 



V. CORRELATION EFFECTS 

In order to compare the numerical simulation results 
with the hydrodynamic approach we scaled them by the 
charge number Z. In Fig. |](a), the energy per particle 
scaled by Z is shown as a function of the relative number 
of electrons n = N/Z. 

The hydrodynamic result (solid curve) was obtained 
by numerically integrating the scaled chemical potential 
(Ell) over the number of electrons n, namely, 



E(r 



ZN ZN 



—— I [idn. 



(26) 



In Fig. ||(b) the same comparison is shown for the scaled 
chemical potential. The numerical simulation results 
were calculated as the difference of the total energies 
for adjacent configurations, namely, /ijv = En+i — Em- 
We see that when the charge number Z increases the 
curves tend towards the hydrodynamic result. This is 
more clearly seen for the chemical potential (Fig. |§|(b)) 
than for the total energy curves (Fig. ||(a)). 

The deviation of the energy and the chemical potential 
of the 'exact' numerical simulation results from their hy- 



drodynamic counterparts has to be interpreted as a corre- 
lation energy, because such correlations are not included 
in the hydrodynamic approach. As an example, we show 
in Fig. |^, for Z = 100, this difference in energy per par- 
ticle and similar results for the chemical potential. The 
curves for other values of the confinement charge demon- 
strate a similar behavior. The most remarkable feature 
of these curves is the linear behavior in the nearly full 
screening region [n — * 1). The curves in the above region 
can be fitted by the rather simple analytical expressions: 



AE/NZ = A E - B E n, 
Afi/Z = — B^n. 

For the results of Fig. || we found 



AE/NZ -- 
Afi/Z 



0.068- 0.034n, 
= 0.065 - 0.067n. 



(27a) 
(27b) 



(28a) 
(28b) 



Although the coefficients of these expressions are in- 
dependent, we found that the ratio of them are Z- 
independent (A E /B E sa 2 and A^/B^ w 1). The lin- 
ear behavior of these curves in the double logarithmic 
plot, as it is seen in Fig. Ilfj, enables us to fit them by 
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FIG. 10: Logarithmic plot of the parameter A of the linear 
fit of the deviation curves for the energy per particle and the 
chemical potential in the nearly complete screening limit. 

A = aZ^ 13 . The straight lines in Fig. [l(] correspond to 
a = 0.69 and (3 = 0.5 for A E , and a = 0.55 and j3 = 0.46 
for Afj,. Therefore, we suggest the following asymptotic 
behaviour for N — > Z 



AE/N ~ \/Z ( 1 



Afj, 



N 
2Z 



(29a) 
(29b) 



These results express the contribution of correlation to 
the energy and the chemical potential of our Coulomb 
bound classical dot. 

Qualitatively, such dependence can be explained by 
means of the crystallization energy which can be esti- 
mated as 



Err = 



d 2 rp{r)E I {r) 



(30) 



where Ei(r) is the local density crystallization energy 
which we approximated by the result from a homoge- 
neous Wigner crystal |ll| 



(31) 



The coefficient equals ao w 3.921 for the triangular lat- 
tice. Now inserting the above expression into Eq. (30), 
replacing the electron density by its asymptotic expres- 
sion (fj"l|), and taking Eq. (25b) into account we obtained 
the following estimate: 
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(32) 



Taking only the first term into account and dividing the 
above crystallization energy by N = Z{1 — (1 — n)} we ob- 
tained the following approximate crystallization energy 
per particle 
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(33) 
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FIG. 11: Comparison of maximal (R) and mean radii (r) of 
the electron system obtained from the numerical simulation 
(for different Z-values) and in the hydrodynamic approach as 
function of n = N/z. The Z-dependence of the maximum and 
the mean radius is enlarged in the inset. 



which demonstrates t he sa me parameter dependencies as 
obtained earlier, Eq. (29a), from a fitting of our simula- 
tion results. 



Differentiating the above crystallization energy expres- 
sion by n we also obtained an estimation for the correla- 
tion contribution to the chemical potential: 



dN 



: E r 



2a^\fZ d 
hxPh: dn 



n{2 



Z(l - n) (34) 



which coincides with our earlier result, Eq. (29b) 



Unfortunately, no similar simple expressions can be ob- 
tained for the dot radius as depicted in Fig. |ll|. We see 
that the radii from the hydrodynamic approach and the 
discrete system differs substantially, and the deviation 
grows in the limiting n — ► 1 case. The discrete system 
is more compact as compared with the continuous one, 
which is an indication that the discrete system may lead 
to overcharging R] which cannot be described by a hy- 
drodynamic theory. 
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VI. CONCLUSIONS 

We studied numerically the ground state properties of 
a 2D model system consisting of classical charged parti- 
cles which are Coulomb bound. This system is similar to 
an atomic system. The N electrons moving in a plane are 
confined through a positive remote impurity potential of 
charge Z. The configurations for n <C N/Z are similar 
to those of a parabolic confined dot, but for n ~ 1 they 
are very different, with an inner core consisting of ap- 
proximately a triangular lattice and an outer region with 
particles situated on a ring. 

A hydrodynamic analysis of the problem clearly em- 
phasizes the importance of correlation effects between 
the negatively charges particles. Furthermore, analyti- 



cal expressions for the correlation energy contributions 
to the total energy and the chemical potential were ob- 
tained in the limit of nearly overscreening, e. g. n ~ 1. 
Those analytical results compare favorably well with our 
numerical "exact" simulation results. 
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APPENDIX A: COULOMB POTENTIAL 
EXPANSION COEFFICIENTS 



According to Eq. ( jiq ) the Coulomb potential expan- 
sion coefficients A n depend on the single variable b = 
1/R, and they can be calculated straightforwardly. In- 
deed, multiplying Eq. (|l8| ) by P^m (t) , integrating it over 
r and using the normalization integral for symmetric Leg- 
endre polynomials one obtains the following integral ex- 
pression: 



A n {b) = (4n+l) 



P2n(T)dT 

Vl + b 2 -T 



(Al) 



The most accurate way to calculate the above integral is 
to expand the denominator into a r 2 ™ power series and 
use the analytical expression for 



2r(-fc)r(n + fc + 3/2) 



what enables us to convert the integral flAip into the 
following sum 



Mb) 



(An + 1) 



2 v ^F(l + 6 2 )"+ 1 /2 

r 2 (m + n+l/2)(l + 6 2 )- 



E 



r(m + l)r(m + 2n + 3/2) 



(A3) 



The convergence of this sum is rather slow. Fortunately, 
it can be improved by adding the asymptotic, namely, 
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replacing Eq. (A3) by the following expression 
{An + 1) 



A n (b) = 



2 v ^F(l + 6 2 )»+ 1 /2 

r 2 (m + n+l/2)(l + 6 2 )" m 
r(m + l)r(m + 2n + 3/2) 



E 

m— 



4n+ 1 



(1 + b 2 ) s + n+1 / 2 

|=-6e b2s [l-$(6^)]}- (A4) 



where s is an integer which should be optimized for rapid 
convergence and $(z) is the error function. The above 
expression is very convenient and enables us to obtain 
good accuracy in a fast way if one uses the recurrence 
expressions for the calculations of the gamma functions 
(as for the Legendre polynomials in Eq. (|22] ) as well). 



